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DRAFT VERSION 

Massive galaxies in cosmological simulations: UV-selected sample 

at redshift z = 2 

^ ] Kentaro Nagainine^, Renyue Cen^, Lars Hernquist^, Jeremiah P. Ostriker^'^, 

O \ & Volker Springel^ 

(N 

< 

•^ . We study the properties of galaxies at redshift z = 2 in a A cold dark matter 

^ . (ACDM) universe, using two different types of hydrodynamic simulation meth- 

-^ . ods - Eulerian TVD and smoothed particle hydrodynamics (SPH) - and a spee- 

ch . trophotometric analysis in the Un,G,R filter set. The simulated galaxies aX z = 2 

^ . satisfy the color-selection criteria proposed by Adelberger et al. (2004) and Stei- 

O , del et al. (2004) when we assume Calzetti extinction with E{B — V^) = 0.15. 

Q . We find that the number density of simulated galaxies brighter than R < 25.5 at 

2; = 2 is about 2 x 10^^ h^ Mpc"'^ for E{B — V) = 0.15 in our most representative 
run, roughly one order of magnitude larger than that of Lyman break galaxies at 
P ! z = 3. The most massive galaxies at z = 2 have stellar masses > lO^^M©, and 

c/3 ! their observed-frame G — R colors lie in the range 0.0 < G — R < 1.0. They typi- 

cally have been continuously forming stars with a rate exceeding 30 M© yr^^ over 
a few Gyrs from 2; = 10 to 2; = 2, although the TVD simulation indicates a more 
sporadic star formation history than the SPH simulations. Of order half of their 
p^ ', stellar mass was already assembled by 2; ~ 4. The bluest galaxies with colors 

—0.2 < G—R < 0.0 at 2; = 2 are somewhat less massive, with Mgtar < l0^^h~^MQ, 
and lack a prominent old stellar population. On the other hand, the reddest mas- 
sive galaxies at 2 = 2 with G — R > 1.0 and Mgtar > 10^^h~^MQ completed the 
build-up of their stellar mass by z ~ 3. Interestingly, our study indicates that 
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the majority of the most massive galaxies at z = 2 should be detectable at rest- 
frame ultra-violet wavelengths, contrary to some recent claims made on the basis 
of near-infrared studies of galaxies at the same epoch, provided the median ex- 
tinction is less than E{B — V) < 0.3 as indicated by surveys of Lyman break 
galaxies at z = 3. However, our results also suggest that the fraction of stellar 
mass contained in galaxies that pass the color-selection criteria used by Steidel 
et al. (2004) could be as low as 50% of the total stellar mass in the Universe at 
z = 2. Our simulations imply that the missing stellar mass is contained in fainter 
{R > 25.5) and intrinsically redder galaxies. The bright end of the rest-frame V- 
band luminosity function oi z = 2 galaxies can be characterized by a Schechter 
function with parameters (<!>*, M^, a) = (1.8 x 10~^, —23.4, —1.85), while the 
TVD simulation suggests a flatter faint-end slope of a ~ —1.2. A comparison 
with z = 3 shows that the rest-frame l^-band luminosity function has brightened 
by about 0.5 magnitude from z = 3 to z = 2 without a significant change in 
the shape. Our results do not imply that hierarchical galaxy formation fails to 
account for the massive galaxies at z > 1. 

Subject headings: cosmology: theory — stars: formation — galaxies: formation 
— galaxies: evolution — methods: numerical 



1. Introduction 

A number of recent observational studies have revealed a new population of red, massive 
galaxies at redshift z ~ 2 (e.g. Chen et al. 2003; Daddi et al. 2004; Franx et al. 2003; 
Glazebrook et al. 2004), utilizing near-infrared (IR) wavelengths which are relatively less 
affected by dust extinction. At the same time, a number of studies focused on the assembly 
of stellar mass density at high redshift by comparing observational data and semi-analytic 
models of galaxy formation (e.g. Poll et al. 2003; Fontana et al. 2003; Dickinson et al. 
2003a). These works argued that the hierarchical structure formation theory may have 
difficulty in accounting for sufficient early star formation. These concerns grew with the 
mounting evidence for high redshift galaxy formation including the discovery of Extremely 
Red Objects (EROs) at ^ > 1 (e.g. Elston et al. 1988; McCarthy et al. 1992; Hu & Ridgway 
1994; Cimatti et al. 2003; Small et al. 2002), sub- millimeter galaxies aX z > 2 (e.g. Small 
et al. 1997; Chapman et al. 2003), Lyman break galaxies (LBGs) at 2; ~ 3 (e.g. Steidel et al. 
1999), and galaxies at 2; > 4, detected either by their Lyman-a emission (e.g. Hu et al. 1999; 
Rhoads & Malhotra 2001; Taniguchi et al. 2003; Kodaira et al. 2003; Ouchi et al. 2003a) 
or by their optical to near-IR colors (e.g. Iwata et al. 2003; Ouchi et al. 2003b; Dickinson 
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et al. 2003b; Chen & Marzke 2004). We now face the important question as to whether 
this evidence for high-redshift galaxy formation is consistent with the concordance ACDM 
model. 

The redshift range around 2; ~ 2 is a particularly interesting epoch also for another 
reason. The redshift interval 1.4 < z < 2.5 has long been known as the 'redshift desert' 
(Abraham et al. 2004; Steidel et al. 2004) for galaxy surveys, because there was no large 
volume limited sample of galaxies available in this regime until very recently. This is because 
strong emission hues from H II regions of galaxies, such as [O ll] A3727, [O ill] AA 4959, 5007, 
Ha and H/9, redshift out of optical wavelengths above 9300 A. These spectral features are 
necessary to easily identify galaxy redshifts, allowing ground-based telescopes to benefit from 
the low-night sky background and high atmospheric transmission in the wavelengths range 
4000 - 9000 A. 

However, recently Adelberger et al. (2004) and Steidel et al. (2004) have introduced new 
techniques for exploring the 'redshift desert', making it possible to identify a large number 
of galaxies efficiently with the help of a color selection criteria in the color-color plane of 
Un — G vs. G — R. In this technique, galaxies at z = 2 — 2.5 are located photometrically from 
the mild drop in the f/„ filter owing to the Ly-a forest opacity, and galaxies at z = 1.5 — 2 
are recognized from the lack of a break in their observed-frame optical spectra. The large 
sample of galaxies identified by these authors at z = 2 makes it now possible to study galaxy 
formation and evolution for over 10 Gyrs of cosmic time, from redshift z = 3 to 2; = 0, 
without a significant gap. We note that the epoch around z = 2 is particularly important for 
understanding galaxy evolution because at around this time, the number density of quasi- 
stellar objects (QSOs) peaked (e.g. Schmidt 1968; Schmidt et al. 1995; Fan et al. 2001) and 
the ultra-violet (UV) luminosity density began its decline by about an order of magnitude 
from z ~ 2 to 2; = (e.g. Lilly et al. 1996; Connolly et al. 1997; Sawicki et al. 1997; Treyer 
et al. 1998; Pascarelle et al. 1998; Cowie et al. 1999). 

These recent observational studies of galaxies at 2; = 2, both in the UV and near-IR 
wavelengths, imply a range of novel tests for the hierarchical structure formation theory. In 
this study, we analyze the properties of massive galaxies at 2; = 2 formed in state-of-the-art 
cosmological hydrodynamic simulations of the ACDM model, and we compare them with 
observations. 

In an earlier recent study (Nagamine et al. 2004a), we argued that, based on two different 
types of hydrodynamic simulations (Eulerian TVD and SPH) and the theoretical model of 
Hernquist & Springel (2003) (hereafter H&S model), the predicted cosmic star formation rate 
(SFR) density peaks at 2; > 5, with a higher stellar mass density at ;i; = 3 than suggested by 
current observations (e.g. Brinchmann & Ellis 2000; Cole et al. 2001; Cohen 2002; Dickinson 
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et al. 2003a; Fontana et al. 2003; Rudnick et al. 2003; Glazebrook et al. 2004), in contrast 
to some claims to the contrary (Poli et al. 2003; Fontana et al. 2003). We also compared 
our results with those from the updated semi-analytic models of Somerville et al. (2001), 
Granato et al. (2000), and Menci et al. (2002), and found that our simulations and the H&S 
model predicts an earlier peak of the SFR density and a faster development of stellar mass 
density compared to these semi-analytic models. 

It is then interesting to examine what our simulations predict for the properties of 
massive galaxies at 2; ~ 2. In this paper, we analyze for this purpose the same set of 
hydrodynamic simulations that was used in Nagamine et al. (2004a), with a special focus on 
the most massive galaxies at this epoch. In Section 2, we briefly describe the simulations that 
we use. In Section 3, we summarize our method for computing spectra of simulated galaxies 
both in the rest-frame and the observed-frame. In Section 4, we show the color-color diagrams 
and color-magnitude diagrams of simulated galaxies. In Section 5 we discuss the stellar 
masses and number density of galaxies at z = 2. We then investigate the observed-frame 
i?-band luminosity function and the rest-frame \/-band luminosity function in Section 6, 
followed by an analysis of the star formation histories in Section 7. Finally, we summarize 
and discuss the implications of our work in Section 8. 



2. Simulations 

We will discuss results from two different types of cosmological hydrodynamic simula- 
tions. Both approaches include "standard" physics such as radiative cooling/heating, star 
formation, and supernova (SN) feedback, although the details of the models and the param- 
eter choices vary considerably. 

One set of simulations was performed using an Eulerian approach, which employed a 
particle-mesh method for the gravity and the total variation diminishing (TVD) method 
(Ryu et al. 1993) for the hydrodynamics, both with a fixed mesh. The treatment of the 
radiative cooling and heating is described in Cen (1992) in detail. The structure of the code 
is similar to that of Cen & Ostriker (1992, 1993), but the code has significantly improved over 
the years with additional input physics. It has been used for a variety of studies, including the 
evolution of the intergalactic medium (Cen et al. 1994; Cen & Ostriker 1999a,b), damped 
Lyman-a absorbers (Cen et al. 2003), and galaxy formation (e.g. Cen & Ostriker 2000; 
Nagamine, Fukugita, Cen, & Ostriker 2001a,b; Nagamine 2002). 

Our other simulations were done using the Lagrangian smoothed particle hydrodynamics 
(SPH) technique. We use an updated version of the code GADGET (Springel et al. 2001b), 
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which uses an 'entropy conserving' formulation (Springel & Hernquist 2002) of SPH to miti- 
gate problems with energy/entropy conservation (e.g. Hernquist 1993) and overcooling. The 
code also employed a subresolution multiphase model of the interstellar medium to describe 
self-regulated star formation which includes a phenomenological model for feedback by galac- 
tic winds (Springel & Hernquist 2003a), and the impact of a uniform ionizing radiation field 
(Katz et al. 1996; Dave et al. 1999). This approach has been used to study the evolution 
of the cosmic SFR (Springel & Hernquist 2003b), damped Lyman-a absorbers (Nagamine, 
Springel, & Hernquist 2004a,b), Lyman-break galaxies (Nagamine, Springel, Hernquist, & 
Machacek 2004b), disk formation (Robertson et al. 2004), emission from the intergalactic 
medium (Furlanetto et al. 2003, 2004a,b,c,d; Zaldarriaga et al. 2004), and the detectability 
of high redshift galaxies (Barton et al. 2004; Furlanetto et al. 2004e). 

In both codes, at each time-step, some fraction of gas is converted into star particles in 
the regions that satisfy a set of star formation criteria (e.g. significantly overdense, Jeans 
unstable, fast cooling, converging gas flow). Upon their creation, the star particles are tagged 
with physical parameters such as mass, formation time, and metallicity. After forming, they 
interact with dark matter and gas only gravitationally as collisionless particles. 

More specifically, in the TVD simulation, the star formation rate is formulated as 

dt " "* t. ' ^'' 

where c* is the star formation efficiency and t^, is the star formation time-scale. In the TVD 
N864L22 run, c* = 0.075 and t* = max(tdyn, 10'^ yr) were used, where tdyn = v^37r/(32G'ptot) 
is the local dynamical time-scale owing to gravity. At each time-step, a part of the gas in a cell 
is converted to a star particle according to Equation (1), provided the gas is (a) moderately 
overdense {6tot > 5.5), (b) Jeans unstable (mgas > 'mj), (c) cooling fast (tcooi < ^dyn), and 
(d) the ffow is converging into the cell (V • v < 0). 

On the other hand, the SPH simulations analyzed in this study parameterized the star 
formation rate as 

%^(l-/3)ff. (2) 

where (3 is the mass fraction of short-lived stars that instantly die as supernovae (taken 
to be f3 = 0.1). The star formation time-scale t* is again taken to be proportional to the 
local dynamical time of the gas: t*(p) = toipth/ pY^'^, where the value of Q = 2.1 Gyr is 
chosen to match the Kennicutt (1998) law. In the model of Springel & Hernquist (2003a), 
this parameter simultaneously determines a threshold density pth) above which a multiphase 
structure of the gas, and hence star formation, is allowed to develop. This physical density is 
8.6 X lO^/i^M© kpc~ for the simulations in this study, corresponding to a comoving baryonic 
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overdensity of 7.7 x 10^ at z = 0. (See Springel & Hernquist 2003a, for a description of 
how pth is determined self-consistently within the modeL) Note that the density p^ that 
determines the star formation rate in equation (2) is actually the average density of cold 
clouds determined by a simplified equilibrium model of a multiphase interstellar medium, 
rather than the total gas density. However, in star-forming regions, the density is so high 
that the mean density in cold gas always dominates, being close to the total gas density. 

We will further discuss the differences between our star formation recipes in the two 
codes in some more detail in Section 8. The cosmological parameters adopted in the simu- 
lations are intended to be consistent with recent observational determinations (e.g. Spergel 
et al. 2003), as summarized in Table 1 where we list the most important numerical parame- 
ters of our primary runs. While there are many similarities in the physical treatment between 
the two approaches (mesh TVD and SPH), we see that the TVD simulation has somewhat 
higher mass resolution, and the SPH somewhat higher spatial resolution. In that sense, the 
two approaches are complementary and results found in common are expected to be robust. 



3. Analysis Method 

We begin our analysis by identifying galaxies in the simulations as groups of star parti- 
cles. Since the two simulations use inherently different methods for solving hydrodynamics, 
we employ slightly different methods for locating galaxies in the two types of simulations, 
each developed to suit the particular simulation method well. Because the stellar groups of 
individual galaxies are typically isolated, there is little ambiguity in their identification, such 
that systematic differences owing to group finding methods are largely negligible. 

For the TVD simulation, we use a version of the algorithm HOP (Eisenstein & Hut 
1998) to identify groups of star particles. This grouping method consists of two discrete 
steps: in the first step, the code computes the overdensity at the location of each particle, 
and in a second step, it then merges identified density peaks based on a set of overdensity 
threshold parameters specified by the user. There are three threshold parameters: a peak 
overdensity (5peak above which the density peak is identified as a candidate galaxy, an outer 
overdensity Sont which defines the outer contour of groups, and finally a saddle point over- 
density 6sa.d which determines whether two density peaks within the outer density contour 
should be merged. Eisenstein & Hut (1998) suggest using (5out, 4ad, ^pcak) = (80,200,240) 
for identifying dark matter halos in N-body simulations. 

We applied the algorithm to the star particles in our simulation using various combi- 
nations of these parameters, and found that the results do not depend very much on the 
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detailed choice. This reflects the higher concentration and isolation of star particle groups 
as compared to the dark matter, making group finding typically straightforward. However, 
the stellar peak overdensities relative to the mean are so high that the saddle-overdensity 
parameter in HOP is unimportant, and the small galaxies associated with a large central 
galaxy at the center of a large dark matter halo are often not well separated by this algo- 
rithm. Therefore, we decided not to apply the second step of the algorithm. Instead, we 
merge two density peaks identified by the first step only if they are within two neighboring 
cells. By visually inspecting the particle distribution of very massive dark matter halos, we 
found that this method succeeds in separating some of the small galaxies embedded in a large 
dark matter halo especially at the outskirts of the halo, and helps to somewhat mitigate the 
'overmerging problem' in high density regions. However, even with this revised method, the 
central massive galaxy remains afterwards, and the addition of a small number of galaxies 
to the faint-end of the luminosity function does not change the qualitative feature of our 
results. Clearly more work is needed to improve this situation on the grouping in the future, 
for example, using an improved grouping method such as VOBOZ (Neyrinck, Gnedin, & 
Hamilton 2004). 

For the SPH simulations, we employ the same method as in Nagamine et al. (2004b), 
where we identified isolated groups of star particles using a simplified variant of the SUBFIND 
algorithm proposed by Springel et al. (2001a). In short, we first compute an adaptively 
smoothed baryonic density field for all star and gas particles. This allows us to robustly 
identify centers of individual galaxies as isolated density peaks. We then find the full extent 
of these galaxies by processing the gas and star particles in order of declining density, adding 
particles one by one to the galaxies. If all of the 32 nearest neighbors of a particle have 
lower density, this particle is considered to be a new galaxy 'seed'. Otherwise, the particle 
is attached to the galaxy that its nearest denser neighbor already belongs to. If the two 
nearest denser neighbors belong to different galaxies, and one of these galaxies has less than 
32 particles, these galaxies are merged. If the two nearest denser neighbors belong to different 
galaxies and both of these galaxies have more than 32 particles, then the particle is attached 
to the larger group of the two, leaving the other one intact. Finally, we restrict the set of 
gas particles processed in this way to those particles which have at least a density of 0.01 
times the threshold density for star formation, i.e. pth = 8-6 x lO^/i^M0kpc~^ (see Springel 
& Hernquist 2003a, for details on how the star formation threshold is determined). Note 
that we are here not interested in the gaseous components of galaxies - we only include gas 
particles because they make the method more robust owing to the improved sampling of the 
baryonic density field. 

We found that the above method robustly links star particles that belong to the same 
isolated galaxy in SPH simulations. A simpler friends-of-friends algorithm with a small 



linking length can achieve a similar result, but the particular choice one needs to make for 
the linking length in this method represents a problematic compromise, either leading to 
artificial merging of galaxies if it is too large, or to loss of star particles that went astray 
from the dense galactic core, if it is too small. Note that, unlike in the detection of dark 
matter substructures, no gravitational unbinding algorithm is needed to define the groups of 
stars that make up the galaxies formed in the simulations. We consider only galaxies with 
at least 32 particles (star and gas particles combined) in our subsequent analysis. 

For both simulations, each star particle is tagged with its mass, formation time, and 
metallicity of the gas particle that it was spawned from. Based on these three tags, we 
compute the emission from each star particle, and co-add the fiux from all particles for 
a given galaxy to obtain the spectrum of the simulated galaxy. We use the population 
synthesis model by Bruzual & Chariot (2003) which assumes a Chabrier (2003b) initial mass 
function (IMF) within a mass range of [0.1, 100] Mq, as recommended by Bruzual & Chariot 
(2003). Spectral properties obtained with this IMF are very similar to those obtained using 
the Kroupa (2001) IMF, but the Chabrier (2003b) IMF provides a better fit to counts 
of low- mass stars and brown dwarfs in the Galactic disk (Chabrier 2003a). We use the 
high resolution version of the spectral library of Bruzual & Chariot (2003) which contains 
221 spectra describing the spectral evolution of a simple stellar population from t = to 
t = 20 Gyr over 6900 wavelength points in the range of 91 A- 160 /xm. 

Once the intrinsic spectrum is computed, we apply the Calzetti et al. (2000) extinction 
law with three different values of E{B — V)= 0.0, 0.15, 0.3 in order to investigate the impact 
of extinction. These values span the range of extinction estimated from observations of LBGs 
at z = 3 (Adelberger & Steidel 2000; Shapley et al. 2001). Using the spectra computed in this 
manner, we then derive rest-frame colors and luminosity functions of the simulated galaxies. 

To obtain the spectra in the observed frame, we redshift the spectra and apply absorption 
by the intergalactic medium (IGM) following the prescription of Madau (1995). Once the 
redshifted spectra in the observed frame are obtained, we convolve them with different filter 
functions, including f/„, G, R (Steidel & Hamilton 1993) and standard Johnson bands, and 
compute the magnitudes in both AB and Vega systems. Apparent f/„, G, R magnitudes are 
computed in the AB system to compare our results with Adelberger et al. (2004) and Steidel 
et al. (2004), while the rest-frame V-band magnitude is computed in the Vega system. 
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4. Color-color & color-magnitude diagrams 

In Figure 1, we show the color-color diagrams of simulated galaxies at z = 2 in the 
observed-frame Un — G vs. G — R plane, for both SPH (top 2 panels) and TVD (bottom 
right panel) simulations. Only galaxies brighter than R = 25.5 are shown in order to match 
the magnitude limit of the sample of Steidel et al. (2004). We plot different symbols for three 
different values of Calzetti extinction: E{B — V)= 0.0 (blue dots), 0.15 (green crosses), and 
0.3 (red open squares). As the level of extinction by dust is increased from E{B — V) = 0.0 to 
0.3, the measured points move towards the upper right corner of each panel. This behavior 
is expected for a conventional star-forming galaxy spectrum, as demonstrated in Figure 2 of 
Steidel et al. (2003). The color-selection criteria used by Adelberger et al. (2004) and Steidel 
et al. (2004) is shown by the long-dashed boxes. The upper (lower) box corresponds to their 
BX (BM) criteria for selecting galaxies at z = 2.0 — 2.5 {z = 1.5 — 2.0). The number density 
of all the galaxies shown (i.e. R < 25.5) is given in each panel for E{B — V) = 0.0, 0.15, 
0.30 from top to bottom, respectively. 

It is encouraging to see that in all the panels most of the simulated galaxies actually 
satisfy the observational color selection criteria for the case of E[B — V) = 0.15. This 
suggests that the simulated galaxies have realistic UV colors compared to the observed ones. 
Simulated galaxies with no extinction tend to be too blue compared to the color-selection 
criteria. In the SPH G6 run, the distribution is wider than in the SPH D5 and TVD 
N864L22 runs, owing to the larger number of galaxies in a larger simulation box, and the 
distribution actually extends beyond the color selection boundaries. The spatial number 
density of galaxies with R < 25.5 is about 2 x 10"^ h^ Mpc"^ for the SPH G6 run. We will 
discuss the number density of galaxies with various cuts in Figure 4 in more detail. 

We repeated the same calculation with the Sloan Digital Sky Survey's u, g, r, i, z filter 
set, and found that the distribution of points hardly changes, even when we exchange Un — G 
with u — g, and G — R with g — r. This means that searches for galaxies in the redshift range 
1.4 < z < 2.5 can also substitute the Sloan filters u,g,r when applying a color-selection 
following the 'BX' and 'BM' criteria adopted by Adelberger et al. (2004) and Steidel et al. 
(2004). 

In Figure 2, we show the simulated galaxies at z = 2 in the plane of apparent R 
magnitude and G — R color. Again, we use three different symbols for three different values 
of extinction: E{B — V)= 0.0 (blue dots), 0.15 (green crosses), and 0.3 (red open squares). 
The vertical long-dashed line and the arrow indicate the magnitude limit oi R = 25.5 used 
by Steidel et al. (2004). The horizontal line at G — i? = 0.75 corresponds to the upper limit 
of the color-selection box of Steidel et al. (2004) and Adelberger et al. (2004). 
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We see that most of the galaxies brighter than R = 25.5 automatically satisfy the 
criterion G — R < 0.75, and almost no galaxies with R < 25.5 fall out of the region for 
E{B — V) = 0.15. There is a significant population of dim {R > 27) galaxies with G — R > 1. 
We will see below that these are low-mass galaxies with stellar masses Mgtar < IO^^H^^Mq. 
We will discuss the amount of stellar mass contained in the galaxies that satisfy the color 
selection criteria in Section 5. 



5. Galaxy stellar masses at z = 2 

In Figure 3, we show the apparent Rab magnitude vs. stellar mass of simulated galaxies 
at redshift 2; = 2, for both SPH (top 2 panels) and TVD (bottom right panel) simulations. 
The three different symbols correspond to three different values of extinction: E[B — V) = 0.0 
(blue dots), 0.15 (green crosses), and 0.3 (red open squares). The vertical long-dashed line 
and the arrow indicate the magnitude limit of R = 25.5 used by Steidel et al. (2004). 

From this figure, we see that there are many galaxies with stellar masses larger than 
10^^h~^MQ and R < 25.5 at 2; = 2 in our simulations. In the SPH G6 run, the masses of 
the most luminous galaxies are substantially larger than those in the D5 run. But this is 
primarily a result of a finite box size effect; the brightest galaxies have a very low space- 
density, and so they are simply not found in simulations with too small a volume. In the 
SPH G6 run, the most massive galaxies at z = 2 have masses 10^"*^ < Mgtar < 10^'^h^^MQ 
with a number density of 3.5 x lO^^/i^ Mpc~^. Such galaxies are hence not commonly found 
in a comoving box-size of ~ (30Mpc)^, in fact, only one such galaxy exists in the D5 run. 
Cosmic variance hence significantly affects the resulting number density, and one ideally 
needs a simulation box larger than ~ (100/i~^ Mpc)^ to obtain a reliable estimate of the 
number density of such massive galaxies. We note however that a couple of such massive 
galaxies are also found in the TVD N864L22 run, which is presumably in part a result of 
overmerging owing to the relatively large cell size compared to the gravitational softening 
length of the D5 run. 

In Figure 4, we summarize the results on the cumulative number densities of galaxies, 
either measured above a certain threshold value of stellar mass [panel (a)], or below a thresh- 
old value of observed-frame R-hand magnitude [panel (b)], or rest-frame V-band magnitude 
[panel (c)]. In panel (a), we see that the results of different runs agree reasonably well in 
the threshold mass range of 10^ < Mstar < lO^^/i^^M©, but differ somewhat at the lower and 
higher mass thresholds. At the low-mass end, the lower resolution of the G6 run results in a 
lower number density than in the D5 run. At the high-mass end, the smaller box-size limits 
the number density in the D5 run compared to the higher value in the 06 run. The result of 
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the TVD run falls in between those of D5 and G6 at the low-mass end, but gives larger values 
relative to the SPH runs at the high-mass end, probably owing to the overmerging problem. 
A similar level of agreement can be seen in panel (c), where the rest-frame magnitude is used 
for the a;-axis, but in panel (b), comparatively large differences are present partly owing to 
the larger stretch in the abscissa. 

We find that the amount of stellar mass contained in the galaxies that satisfy the 
color-selection criteria is in general far from being close to the total stellar mass. Owing 
to the limited box-size and resolution effects, it is somewhat difficult to obtain an accurate 
estimate of the mass fraction above some magnitude limit. However, we list the stellar 
mass fractions as a reference for future work. In the G6 run, which has comparatively low 
mass resolution, the stellar mass fraction contained in galaxies that pass the criteria (i.e. 
R < 25.5 & color-selection) account for only (60, 65, 42)% of all the stars, for extinctions 
of E{B — V) = (0.0,0.15,0.3), respectively. These fractions are relatively large compared 
to other runs, because, as we will see in Section 6, the luminosity function in the G6 run 
starts to fall short of galaxies near R ~ 25.5 compared to the D5 run owing to its limited 
resolution, and the G6 run lacks low-mass faint galaxies with R > 25.5. On the other hand, 
in the D5 run, the stellar mass fraction contained in galaxies that pass the criteria (i.e. 
R < 25.5 & color-selection) account for only (32, 17, 2)%, and (16, 13, 1)% of the stars in 
the TVD N864L22 run, for extinctions of E{B - V) = (0.0,0.15,0.3), respectively. These 
fractions are smaller compared to those of the G6 run, because the D5 and TVD N864L22 
runs lack massive bright galaxies owing to their smaller box sizes. If we relax the criteria 
from 'i? < 25.5 & color-selection' to 'only R < 25.5', then the corresponding stellar mass 
fractions increase to (75, 65, 44)% for G6, (32, 17, 3)% for D5, and (46, 34, 29)% for TVD 
N864L22. In summary, we expect that the true stellar mass fraction that satisfies the criteria 
lies somewhere between the results of G6 and D5 (and TVD N864L22) runs, which should 
be roughly ~ 50%. 

This implies that the current surveys with a magnitude limit of i? ~ 25.5 at best 
account for half of the stellar mass in the Universe. Nagamine et al. (2004a) reached a 
similar conclusion based on different theoretical arguments that compared the results of our 
hydrodynamic simulations and the theoretical model of Hernquist & Springel (2003) with 
near-IR observations of galaxies (e.g. Cole et al. 2001; Dickinson et al. 2003a; Fontana et al. 
2003; Rudnick et al. 2003; Glazebrook et al. 2004). Franx et al. (2003) and Daddi et al. 
(2004) suggested that the unaccounted-for contribution to the stellar mass might be hidden 
in a red population of galaxies. Our simulation suggests that the missed stellar masses are 
mostly contained in the fainter {R > 25.5) and redder galaxies. But note that we have not 
considered large values of extinction with E{B — V) > 0.3. 
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In Figure 5, we show the observed- frame G — R color vs. stellar mass of simulated 
galaxies at z = 2 for both SPH (top 2 panels) and TVD (bottom right panel) simulations. 
Here, only the case of E{B — V)= 0.15 is shown. In the top two panels for the SPH runs, 
the red open crosses show the galaxies that are brighter than R = 25.5, and the blue dots 
show the rest of the galaxies. The discrete stripes seen at low-mass end of the distribution 
for the SPH runs are due to the discreteness of the star particles in the simulation. For the 
TVD run, there is a strong concentration of points at G — i? ~ 1.4, which we will discuss in 
detail in Section 7. 

An interesting point to note in this figure is that the bluest galaxies aX z = 2 are not 
necessarily the most massive ones. Rather, the bluest galaxies with colors —0.2 < G—R < 0.2 
are galaxies with somewhat lower mass in the range Mstar < 10^'^/i~^Mq. The most massive 
galaxies with Mstar ^ IO^^/i^^Mq have slightly redder colors of 0.2 < G — R < 0.6 owing to 
their underlying old stellar population. 

In the TVD simulation, there are a couple of galaxies with G — R > 2.0 and Mstar > 
10^°/i^^Mq, which are absent in the SPH simulations. These galaxies are dominated by the 
old stellar population because they have completed their star formation before z ~ 3 as we 
will see in Section 7. 



6. Galcixy luminosity functions 

Figure 6 shows the observed-frame R-hand luminosity function for all simulated galaxies 
at z = 2 in both SPH (top 2 panels) and the TVD (bottom left panel) simulations. In 
these 3 panels, blue long-dashed, green solid, and red dot-short-dashed lines correspond to 
E{B — V) = 0.0, 0.15, 0.30, respectively. The vertical dotted lines indicate the magnitude 
limit oi R = 25.5 used by Steidel et al. (2004). In the bottom right panel, a combined 
result is shown with the SPH D5 (green solid line), G6 (red long-dashed line), and TVD 
(blue dot-dashed line) results overplotted for the case of E{B — V)= 0.15. For comparison, 
Schechter functions with the following parameters (simply chosen by eyeball fit) are shown 
as short-dashed lines: (<!>*, M^, a) = (1 x 10-^,23.2, -1.4) for SPH D5 & G6 runs (top two 
panels) and the combined results (bottom right panel), and (1 x 10^^,24.5, —1.4) for TVD 
(bottom left), respectively. 

The bottom right panel of Figure 6 gives a relative comparison between the different 
runs. We see that the SPH G6 run contains a larger number of luminous galaxies with 
i? < 25 compared to the other runs, an effect that results from its larger box-size. On 
the other hand, the G6 run lacks fainter galaxies at i? > 25 compared to D5 run owing to 
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its limited mass resolution. The SPH G6 run hence covers the bright-end of the Schechter 
function with M^ = 23.2, while the D5 run covers better the faint-end of the luminosity 
function with a slope of a = —1.4. The TVD N864L22 run agrees with the D5 results at the 
bright-end near R = 25, but the deviation from the D5 result starts to become somewhat 
large at i? > 28. We will discuss a number of possible reasons for these differences between 
the SPH and TVD runs in Section 8. 

In Figure 7, we show the rest-frame l^-band luminosity function of all simulated galaxies 
at z = 2 for both SPH (top 2 panels) and TVD (bottom left panel) simulations. In these 
three panels, blue long-dashed, green solid, and red dot-short-dashed lines correspond to 
extinctions of E{B — V") = 0.0, 0.15, 0.30, respectively. In the bottom right panel, the 
combined result is shown with the SPH D5 (green solid line), G6 (red long-dashed), and 
TVD (blue dot-dashed) runs overplotted for the case of E{B — V) = 0.15. For reference, 
the magnitude limit of My = —21.2 (for h = 0.7) for the survey oi z = 3 LBGs by Shapley 
et al. (2001) is indicated by the vertical dotted line and the arrow. For comparison, we 
also include Schechter function fits with parameters ($*, My, a) = (1.8 x 10^^, —23.4, —1.8) 
as short-dashed curves in all panels. These Schechter parameters are the same as the ones 
that Shapley et al. (2001) fitted to their observational data, except that the characteristic 
magnitude My is 0.5 magnitudes brighter. The same amount of brightening is seen in the 
G6 run from z = 3 (see Fig. 5 of Nagamine et al. (2004b) ^) to z = 2. So the G6 run with 
E{B — V) = 0.15 describes the bright-end of the luminosity function aX z = 2 and z = 3 
quite well. In the bottom left panel for the TVD run, another Schechter function with a 
much shallower faint-end slope and parameters {^*,My,a) = (3.5 x 10^^,-22.5,-1.15) is 
also shown with a magenta long-dash-short-dashed line for comparison. One can see that 
the TVD run contains more fainter galaxies with My > —15 than the SPH runs. This could 
be due to the higher baryonic mass resolution and the additional small-scale power in the 
TVD run compared with the SPH runs as we will discuss more in Section 7. 

We can compute the expected characteristic magnitude M^ oi z = 2 galaxies from the 
observational estimate of M|j at z = 3. The difference in the luminosity distance between 
z = 3 and z = 2 corresponds to about one magnitude for our adopted cosmology, and we 
saw in the previous paragraph that the rest-frame \^-band luminosity function brightened 
by ~ 0.5 magnitude from z = 3 to z = 2. Therefore, we expect that M^ should be brighter 
by about 1.5 magnitudes at z = 2 compared to z = 3. Adelberger & Steidel (2000) reported 



^We note that there was an error in the calculation of the AB magnitudes in Nagamine ct al. (2004b). 
An additional factor of 51og(l + z) = 3.01 has to be added to obtain correct AB magnitudes. The corrected 
version of the paper will be posted at astro-ph/0311295. The rest-frame y-band results were not affected 
by this error. 
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M^ = 24.54 for the observed-frame i?-band luminosity function of the z = 3 galaxies. 
Therefore we obtain M^ ~ 23.0 as the expected value for the z = 2 galaxies. The result of 
the SPH G6 run shown in Fig. 6 is reasonably close to this expected value. It is also clear 
that the box size of the D5 and TVD runs are too small to reliably sample the brightest 
galaxies at z = 2. This is also seen in Fig. 7 as a lack of brightest galaxies at My < —23 in 
the D5 and TVD runs. This is because the brightest galaxies form in rare high density peaks 
which are at the knots of the large-scale filaments, and one needs a sufficiently large volume 
to sample a fair number of such bright galaxies reliably. Since the comoving correlation 
length of LBGs aX z = 3 are ~ 4/i~^ Mpc (Adelberger et al. 2003), one needs a box size 
larger than at least Lbox ^ 40/i~^ Mpc in order to have a fair sample of LBGs aX z = 3. 
The simulated volume of the G6 run is (100/i~^ Mpc)^ comoving, and one of the face of the 
simulation box corresponds to 2.5 deg^ at z = 3. The coverage area of the LEG survey of 
Steidel et al. (2003) at z = 3 is about 0.4 deg^, and the depth of their survey is Az ~ 0.5 
which corresponds to ~ 500h~^ Mpc comoving for our adopted flat-A cosmology. Therefore 
the volumes of the G6 simulation and the Steidel et al. survey are comparable. 



7. Star formation history of galaxies 

In Figures 8, 9, and 10, we show the star formation (SF) histories of massive galaxies 
at 2; = 2 as a function of cosmic age, for the SPH D5, G6, and TVD runs, respectively. 
In these figures, panels (a) & (6) show results for the two most massive galaxies in each 
simulation box, while panels (c) & (d) are for the two reddest galaxies among those with 
M^tar > IO^^/i^^Mq. In the right side of each panel, the following quantities for the case of 
E{B — V) = 0.15 are listed from top to bottom: stellar mass in units of /i^^Mq, rest-frame 
V-band magnitude, observed- frame i?-band magnitude, and G — R color. In panel (e), the 
cumulative star formation history of the galaxies in panels (a) — (d) is shown. 

The results of the SPH D5 and G6 runs (Fig. 8 & 9) suggest that the most massive 
galaxies have almost continuously formed stars with a rate exceeding 30 Mq yr^^ over a few 
Gyrs from z = 10 to z = 2 (note that 30 M© yr"^ x 3 Gyr ~ lO^^M©). Half of their stellar 
mass was already assembled hj z = 4. The reddest galaxies (panels (c) & (d)) have colors 
of G — i? > 1.0, and their star formation has become less active after z = 3, with the bulk 
of their stellar mass created before z = 3 as can be seen in the panel (e). 

The TVD simulation (Fig. 10) on the other hand suggests a more sporadic star formation 
history than the SPH simulations, and the two most massive galaxies in the simulation box 
experienced massive star formation between z = 3 and z = 2 that exceeded 1000 M© yr~^ 
for brief periods. The peak star formation rate of the most massive galaxy in the G6 run 
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also reaches 1000 M© yr^^ occasionally. In the real universe, such violent starburst activity 
is found in high-redshift quasars and submm sources where the existence of dense, warm, 
and massive molecular clouds of 10^'^ — lO^^M© is inferred from the observed rotational 
emission line of CO (e.g. Omont et al. 2003; Beelen et al. 2004). The number density of 
quasar host galaxies with such a violent star formation activity is not well-constrained yet, 
but it probably is not very large because the space density of quasars with Mb < —26 is 
smaller than nq ~ lO^^Mpc"^ (e.g. Pei 1995). The space density of submm sources is also 
similar: nsubmm ~ 6.5 x lO^^Mpc^^ aX z ^ 2 (Chapman et al. 2003). These values are 
much smaller than the space density of Lyman break galaxies nLBG ~ 4 x 10""^ /i^ Mpc~^ 
(Adelberger et al. 2004). Finding 3 such massive starbursts at z = 2 in our TVD simulation 
box of {22h~^ Mpc)^ results in a number density of ~ 3 x lO~^h^M'pc^^, which may be too 
high compared to the observed value. It is also possible that the two galaxies shown in 
Figure 10 are affected by the overmerging problem which tends to merge galaxies in very 
high density regions owing to a lack of spatial resolution. If these two massive galaxies were 
broken up into several less massive galaxies, then the 1000 M© yr~^ of star formation would 
be distributed to a few hundreds Mq yr"-*^ for each galaxy. 

The two reddest massive galaxies in the TVD run (Fig. 10c & lOd) have cumulative star 
formation histories which are similar to those of the SPH runs, with half of their stellar mass 
already assembled by z ~ 4. They have completed their star formation by z ~ 3 and have 
been quiet from z = 3 to z = 2. Because they are dominated by an old stellar population, 
they have very red colors oi G — R> 2.0 and show up in the upper right corner of the TVD 
(bottom right) panel in Fig. 1. 

Another way of looking at the star formation history of galaxies is to look at the dis- 
tribution of formation times of stellar particles in the simulation for a given set of galaxies 
as a whole; i.e. we here analyze the combined star formation history for a particular class 
of galaxies divided by their stellar mass. In Figure 11, we show the mass- weighted proba- 
bility distribution of stellar masses as a function of formation time (tform = corresponds 
to the Big Bang), for different samples of galaxies. To this end, we have split the total 
galaxy sample into 3 categories by stellar mass ai z = 2: Msmr > IO^^H'^Mq (top panel), 
10^ < Mstar < lO^o/i-^M© (middle panel), and Mstar < lO^h'^MQ (bottom panel). Differ- 
ent line types correspond to different simulations: red long-dashed (SPH G6 run), blue solid 
(SPH D5 run), and black dot-dashed (TVD N864L22 run). The vertical dotted lines indicate 
four different epochs, namely z = 2, 3, 4, & 5, as indicated in the top panel. 

In the two top panels for massive galaxies, the distribution roughly agrees between the 
different runs, except for the large bump in the TVD run close to z = 2, which is caused 
by starbursts in a few massive galaxies. Also note that the distribution of G6 run is slightly 
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more bumpy owing to lack of mass resolution compared to other runs. 

In the bottom panel for the less massive galaxies, the distribution is very different 
between the 3 runs. In the TVD run, the bulk of stars in low-mass galaxies forms before 
z = 4, whereas the opposite is true for the SPH G6 run. These low-mass galaxies that formed 
very early on appeared as a strong concentration of points at G — i? ~ 1.2 in Fig. 5 for the 
TVD N864L22 run. The result for the SPH D5 appears to be intermediate between the two. 
This could perhaps be caused by the additional small-scale power in the TVD run compared 
with the SPH runs. The initial mean inter-particle separations of both dark matter and gas 
particles in the two SPH runs are comoving 104 and 206 h~^ kpc for the D5 and G6 run, 
respectively, but the mean inter-particle separation of dark matter particles in the TVD run 
is comoving 51 h'^ kpc, with a baryonic cell size of 25 h^^ kpc. Therefore, at very early times, 
the TVD run should be able to resolve smaller matter fluctuations than the two SPH runs. 
Such an effect could explain the earlier formation epoch of low-mass galaxies in the TVD 
run compared to the two SPH runs, provided it is not compensated by the comparatively 
lower gravitational force resolution (roughly 2 cells) of the TVD particle-mesh method. 



8. Discussion & Conclusions 

We have used two different types of hydrodynamic cosmological simulations (Eulerian 
TVD and SPH) to study the properties of massive galaxies at 2; = 2 in a ACDM universe, 
with particular emphasis on an observationally inspired selection based on the Un,G,R filter 
set. The simulated galaxies aX z = 2 satisfy well the color-selection criteria proposed by Adel- 
berger et al. (2004) and Steidel et al. (2004) when we assume Calzetti et al. (2000) extinction 
with E{B — V) = 0.15. However, we find that the fraction of stellar mass contained in galax- 
ies that pass the color-selection criteria could be as low as 50% of the total. The number 
density of simulated galaxies brighter than R < 25.5 at 2; = 2 is about 2 x 10~^ h^ Mpc~^ 
for E{B — V) = 0.15 in the most representative run (SPH G6 run), roughly one order of 
magnitude larger than that of Lyman break galaxies aX z = 3. The increase of the number 
density of bright galaxies can be viewed as a consequence of the ongoing hierarchical build-up 
of larger and brighter galaxies with ever more stellar mass. 

The most massive galaxies aX z = 2 have stellar masses of > lO^'^^M©, and their observed- 
frame G — R colors lie in the range 0.0 < G — R < 1.0 provided that the median extinction 
is E{B — V) ~ 0.15. They have been continuously forming stars with a rate exceeding 
30 M0 yr~^ over a few Gyrs from z = 10 to z = 2. Typically, half of their stellar mass 
was already assembled by z ~ 4. The bluest galaxies with —0.2 < G — R < 0.0 at z = 2 
are somewhat less massive, with Mgtar < l0^^h~^MQ, and are less dominated by old stellar 
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populations. Our study suggests that the majority of the most massive galaxies aX z = 2 
could be detected at rest-frame UV wavelengths contrary to some recent claims based on 
near-IR studies (e.g. Franx et al. 2003; van Dokkum et al. 2004) of galaxies at the same 
epoch, provided the median extinction is in the range E{B — V) < 0.3 as indicated by the 
surveys of Lyman break galaxies at ^ = 3. We plan to extend our analysis to the near-IR 
pass-bands in subsequent work. We find that these massive galaxies reside in the highest 
overdensity regions in the simulations at z = 2. 

While it is encouraging to see that the two very different types of hydrodynamic simu- 
lation methods give a similar picture for the overall properties of massive galaxies aX z = 2, 
we have also seen some notable differences in the results of the SPH and TVD simulations 
throughout the paper. 

For example, the star formation history of galaxies in the TVD simulation is clearly 
more sporadic or 'bursty' than that in the SPH simulations. This probably owes to the 
differences in resolution and the details of how star formation is modeled in the two codes, 
as we described in Section 2. While the dependence of the star formation rate on local gas 
density is essentially the same for the TVD and SPH simulations, the different numerical 
coefficients effectively multiplying the local density on the right hand sides of Equations 
(1) and (2) produce different temporal smoothing of the local star formation rates and the 
parameters adopted by the TVD code leads to a more "bursty" output (pj^°/p^^^ — 3.8), 
although not to a significantly different total star formation rate. At this point we simply 
must await a more detailed theory or be guided by observations in choosing between the two 
approaches. 

Another difference we saw is that the faint-end slope of the rest-frame l^-band luminosity 
function: the TVD simulation gives a = —1.15, similar to what is observed locally a = —1.2 
(e.g. Blanton et al. 2001), but the SPH D5 and G6 runs give much steeper slope of a = —1.8. 
The question of the nature of the faint-end slope of the luminosity function is one that we 
have not resolved at this point. The difference may owe to the difference in the treatment of 
feedback, but Chiu et al. (2001) found a flat faint-end slope oi a = —1.2 in a high-resolution 
softened Lagrangian hydrodynamic (SLH) simulation (Gnedin 1995) even when the code did 
not include SN feedback, suggesting that the photoionization of the gas in low-mass halos by 
the UV background radiation field is responsible for suppressing the formation of low-mass 
galaxies (see also Thoul & Weinberg 1996; Quinn et al. 1996; Bullock et al. 2000). Thus 
the theoretical expectation is not fully clear at present, and it will be very interesting to see 
what is indicated by upcoming observational programs of deeper high redshift galaxies and 
future theoretical studies. 

The present study shows that despite a considerably different numerical methodology 
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and large differences in the assumed model for star formation, two independent hydrody- 
namic codes predict the formation of rather massive galaxies at 2; ~ 2 in the ACDM cosmol- 
ogy, with properties that are broadly consistent with current observations. The results do 
not suggest that hierarchical galaxy formation fails to be able to account for these objects. 
On the contrary, it predicts that significant additional mass in galaxies will be found as the 
observational selection criteria are broadened. 
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Table 1. Simulations 



Run Lbox [h ^ Mpc] A^mcsh/ptci "^dm [h ^Mq] mgas [h ^Mq] M [h ^ kpc] 



TVD: N864L22'^ 


22.0 


8643 


8.9 X 10^ 


2.2 X 10^ 


25.5 


SPH: DS'' 


33.75 


3243 


8.2 X 10^ 


1.3 X 10^ 


4.2 


SPH: G6^ 


100.0 


486^ 


6.3 X 10^ 


9.7 X 10^ 


5.3 



Note. — Parameters of the primary simulations on which this study is based. The quantities 
listed are as follows: Lbox is the simulation box size, A'^tncsh/ptci is the number of the hydrody- 
namic mesh points for TVD or the number of gas particles for SPH, ttidm is the dark matter 
particle mass, mgas is the mass of the baryonic fluid elements in a grid cell for TVD or the 
masses of the gas particles in the SPH simulations. Note that TVD uses 432^ dark matter par- 
ticles for N864 runs. A£ is the size of the resolution element (cell size in TVD and gravitational 
softening length in SPH in comoving coordinates; for proper distances, divide by 1 + z). The 
upper indices on the run names correspond to the following sets of cosmological parameters: 
(fiM,fiA,fib,/i,n,a8) = (0.29,0.71,0.047,0.7,1.0,0.85) for (a), and (0.3,0.7,0.04,0.7,1.0,0.9) 
for (b). 
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Fig. 1. — Color-color diagrams of simulated galaxies at z = 2 in the observed-frame Un — G 
vs. G — R plane for both SPH (top 2 panels) and TVD (bottom right panel) simulations. 
Only those galaxies that are brighter than R = 25.5 are shown to match the magnitude limit 
of the Steidel et al. (2004) sample. The three different symbols represent three different 
values of Calzetti extinction: E{B — V) = 0.0 (blue dots), 0.15 (green crosses), and 0.3 (red 
open squares). The number density of galaxies shown (i.e. R < 25.5) is given in the bottom 
right corner of each panel for E{B — V) = 0.0, 0.15, 0.30 from top to bottom, respectively. 
The color-selection criteria used by Adelberger et al. (2004) and Steidel et al. (2004) are 
shown by the long-dashed boxes. The upper (lower) box corresponds to their BX (BM) 
criteria for selecting out galaxies at z = 2.0 — 2.5 {z = 1.5 — 2.0). 
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Fig. 2. — Color-magnitude diagrams of simulated galaxies at z = 2 on the i? vs. G — R 
plane for both SPH (top 2 panels) and TVD (bottom right panel) simulations. All simulated 
galaxies are shown. Blue dots, green crosses, and red open squares correspond to E{B — V) = 
0.0, 0.15, 0.30. The magnitude limit of i? = 25.5 and the color range of -0.2 <G-R< 0.75 
used by Steidel et al. (2004) and Adelberger et al. (2004) are shown by the long-dashed lines. 
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Fig. 3. — Observed-frame R magnitude vs. stellar masses of simulated galaxies at z = 2 
for both SPH (top 2 panels) and TVD (bottom right panel) simulations. Blue dots, green 
crosses, and red open squares correspond to E{B — V) = 0.0, 0.15, 0.30. The magnitude 
limit oi R = 25.5 used by Steidel et al. (2004) is shown by the vertical long-dashed line. 
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Fig. 4. — Number density of all galaxies with stellar mass [panel (a)] , observed-frame i?-band 
magnitude [panel (b)], rest-frame K-band magnitude [panel (a)] above a certain value on the 
abscissa. The symbols correspond to the results from the SPH D5 run (blue open squares), 
the G6 run (red filled triangles), and the TVD run (green crosses). 
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Fig. 5. — Stellar masses vs. G — R color of simulated galaxies at z = 2 in the observed-frame 
for both SPH (top 2 panels) and TVD (bottom right panel) simulations. The red open 
crosses show the galaxies that are brighter than R = 25.5, while the remainder of the points 
are shown in blue dots. Here, only the case for E{B — V^) = 0.15 is shown. The most massive 
galaxies are not the bluest ones owing to the underlying old stellar population. The color 
range of —0.2 < G — R < 0.75, which is the G — R color range of the selection criteria shown 
in Figure 1, is indicated by the vertical long-dashed lines. 



-30- 







oo 






-3 


CO 




1 
O 


-4 


Oh 




S 




CO 


-5 


^ 





^ 




tJJO 


-1 



o 



3 

■4 

5 




TVD 
N864L22 



Combined 
E(B-V) = 0.1: 




20 



25 



30 



20 



AB 



25 



30 



AB 



Fig. 6. — R band luminosity function in the observed-frame for all simulated galaxies at 2; = 2 
for both SPH (top 2 panels) and TVD (bottom left panel) simulations. For these 3 panels, 
blue long-dashed, green solid, and red dot-short-dashed lines correspond to E{B — V) = 0.0, 
0.15, 0.30, respectively. The vertical dotted line indicates the magnitude limit oi R = 25.5 
used by Steidel et al. (2004). In the bottom right panel, a combined result is shown with 
the SPH D5 (green solid line), G6 (red long-dashed), and TVD (blue dot-dashed) results 
overplotted for the case of E{B — V) = 0.15. For comparison, the Schechter functions with 
following parameters are shown in short-dashed lines: ($*,M|j,, a) = (1 x 10~^,23.2, —1.4) 
for SPH D5 & G6 (top two panels), and combined results (bottom right panel), and (1 x 
10~^,24.5, —1.4) for the TVD run (bottom left panel), respectively. 
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Fig. 7. — Rest-frame l^-band luminosity function of all simulated galaxies at 2 = 2 for both 
SPH (top 2 panels) and TVD (bottom left panel) simulations. For these 3 panels, blue long- 
dashed, green solid, and red dot-short-dashed lines correspond to E{B — V) = 0.0, 0.15, 
0.30, respectively. In the bottom right panel, a combined result is shown with the SPH D5 
(green solid line), G6 (red long-dashed), and TVD (blue dot-dashed) results overplotted for 
the case of E{B - V) = 0.15. The magnitude limit of My = -21.2 (for h = 0.7) for the 
survey of the LBGs at ^ = 3 by Shapley et al. (2001) is shown by the vertical dotted line. To 
guide the eye, a Schechter function with parameters ($*, My, a) = (1.8 x 10~^, —23.4, —1.8) 
is shown in the short-dashed curve in all panels. In the bottom left panel of TVD run, 
another Schechter function with parameters ($*, My, a) = (3.5 x 10^^, —22.5, —1.15) is also 
shown as a magenta long-dash-short-dashed line. 
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Fig. 8. — Panels (a) - (d): Star formation history as a function of cosmic time for selected 
galaxies in the SPH D5 run. In the right side of panels (a) — (d), the following values 
for the case of E{B — V) = 0.15 are listed: stellar mass in units of h'^M^, rest-frame 
V-band magnitude, observed-frame i?-band magnitude, and G — R color. Galaxies (a) & 
(b) are the two most massive galaxies, and (c) & (d) are the two reddest galaxies with 
Mstar > 1 X lO^°/i"-^M0. Panel (e): Cumulative SF history of galaxies (a) - (d). 
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Fig. 9. — Panels (a)-(d): Star formation history as a function of cosmic time for selected 
galaxies in the SPH G6 run. In the right side of each panel, the following values are listed: 
stellar mass in units of /i~^M0, rest-frame V-band magnitude, observed- frame i?-band mag- 
nitude, and G — R color. Galaxies (a) & (b) are the two most massive galaxies, and (c) & 
(d) are the two reddest galaxies with Mstar > 1 x 10^°/i~"'^Mq. Note that panels (a) & (b) 
use a logarithmic scale for the ordinate, while panels (c) & (d) use a linear scale. Panel (e): 
Cumulative SF history of galaxies (a) - (d). 
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Fig. 10. — Panels (a)-(d): Star formation history as a function of cosmic time for selected 
galaxies in the TVD N864L22 run. In each panel, the following values are listed: stellar mass 
in units of /i~^M0, rest-frame V-band magnitude, observed- frame i?-band magnitude, and 
G — R color. Galaxies (a) & (b) are the two most massive galaxies, and (c) & (d) are the two 
reddest galaxies with Mgtar > 1 x 10^^h~^MQ. Note that panels (a) & (b) use a logarithmic 
scale for the ordinate, while panels (c) & (d) use a linear scale. Panel (e): Cumulative SF 
history of galaxies (a) - (d). 
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Fig. 11. — Distribution of stellar masses in the simulations as a function of formation time 
(inform = corresponds to the Big Bang). The galaxy sample is divided into 3 categories 
by their stellar mass: Mstar > IO^^/i^^Mq (top panel), 10^ < Mstar < W°h-^MQ (middle 
panel), Mstar < 10^h~^MQ (bottom panel). Different line types correspond to different 
simulations: red long-dashed (SPH G6 run), blue solid (SPH D5 run), and black dot-dashed 
(TVD N864L22 run). The vertical dotted hues indicate four epochs oi z = 2, 3, 4, & 5, as 
indicated in the top panel. 



